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A new approach to model order reduction of the Navier-Stokes equations at high Reynolds 
number is proposed. Unlike traditional approaches, this method does not rely on empirical 
turbulence modeling or modification of the Navier-Stokes equations. It provides spatial 
basis functions different from the usual proper orthogonal decomposition basis function 
in that, in addition to optimally representing the training data set, the new basis func- 
tions also provide stable and accurate reduced-order models. The proposed approach is 
illustrated with two test cases: two-dimensional flow inside a square lid-driven cavity and 
a two-dimensional mixing layer. 
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1. Introduction 



Numerical simulation of fluid flows is usually a very computationally intensive en- 
deavor. Even when adequate computational resources are available, numerical simula- 
tions often provide too little understanding of the solutions they produce. There are 
significant scientific and engineering benefits in developing and studying Model Order 
Reduction (MOR) techniques capable of producing Reduced Order Models (ROMs) of 
complex fluid flows that retain physical fidelity while substantially reducing the size and 
cost of the computational model. 

In fluid flow applications, Proper Orthogonal Decomposition (POD) and Galerkin pro- 



jection form a popular MOR strategy (Noack et al. 2011 Holmes et al. 2012 1. Despite 



its popularity, deriving robust and accurate ROMs using the POD-Galerkin approach 
remains challenging and many important issues remain unresolved. In particular, appli- 
cation of the POD-Galerkin strategy to turbulent fluid flows remains an active area of 
research. Turbulence is a flow regime characterized by chaotic, multi-scale dynamics in 
both space and time. In the limit of high Re numbers, the dynamics of turbulence are 
qualitatively universal: large scale flow features are broken down into smaller and smaller 



scales until the scales are fine enough that viscous forces can dissipate their energy (Pope 
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2000 Moin & Mahesh 1998 Tennekes & Lumley 1972 ). Application of the standard POD- 



Galerkin strategy to a turbulent flow is problematic because POD, by construction, is 
biased toward the large, energy containing scales of the turbulent flow. ROMs gener- 
ated using only the first most energetic POD basis functions are, therefore, not endowed 
with the natural energy dissipation of the smaller, lower energy turbulent scales. Cur- 
rently, two broad strategies for adapting the POD-Galerkin approach to high Re number 
flows are available: (1) direct modeling of small scale dynamics by inclusion of a suffi- 
ciently large number of POD basis functions and (2) approximate modeling of small scale 
dynamics using empirical turbulence models. Unfortunately, neither strategies is ideal. 
Inclusion of a large number of POD basis functions creates very large POD-Galerkin 
ROMs that are still very computationally expensive to solve. Addition of turbulence 
models is equally undesirable because the empirical terms modify the dynamics of the 
Navier-Stokes equations (|Rempfer fc Fasel||1994| |Aubry et al.||1988| |Ukeiley et aT||2001 



Sirisup fc Karniadakis|2004| |Iliescu fc Wang||2012| |Bailon-Cuba et aZ.||2012| ). 



In this paper, a new approach to MOR of turbulent Navier-Stokes fluid flows is pro- 
posed. Instead of modeling the small, energy-dissipative scales using empirical eddy vis- 
cosity models, the small scales are modeled directly using spatial basis functions that are 
different from the standard POD basis functions. In addition to optimally representing 
the solution, these new basis functions also resolve a greater proportion of the small, 
energy-dissipative scales of the turbulent flow. The proposed methodology is formulated 
as a small-scale constrained minimization problem that can be solved numerically using 
standard, off-the-shelf MATLAB algorithms. 

This paper is organized as follows. In §2 the POD and Galerkin methods are summa- 
rized and the approach is applied to the incompressible Navier-Stokes equations. In §3 
the proposed new approach is introduced and the algorithm and its numerical implemen- 
tation are outlined. In §4 two classical benchmarks, the lid-driven cavity and a mixing 
layer are modeled using the standard POD-Galerkin approach and the new approach. 
Finally, in §5, the main results are summarized and future prospects laid out. 



2. Methodology 

2.1. Spectral methods 

Consider a dynamical system that evolves in a Hilbert space H, u(t) G H governed by 

u = X(u) (2.1) 

where A is a vector field on H. Often, the governing variable u is spatio-temporal, 
u = u(x,t), x € f2, t € [0, T] and thus, A is a spatial differential operator. In the 
spectral approach ( Boyd|200T| Canuto et aL||199ij 2006), the governing variable, u(x,t) 



is discredited using separable basis functions, a,(t) and Ui(x) 

n 

u(x, t) w u [1 -' n] = ai(t)ui(x) (2.2) 

i 

where {it.; € H\i — 1, . . . , n} is a basis for the subspace of H. In the method of lines, the 
spatial basis functions, itj are known a priori and the goal, therefore, is to find temporal 
coefficients, di that satisfy the differential equation. In general, the origin and/or form 
of the spatial basis functions Ui is quite arbitrary. In the context of spectral methods 
in Computational Fluid Dynamics (CFD) the spatial basis functions, Ui are usually 
analytical functions such as trigonometric functions and power functions, i.e. Chebyshev 
polynomials. The benefit of using these functions is that their spatial derivatives are 
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known and numerically efficient algorithms such as the Fast Fourier Transform (FFT) 
can be utilized. In the context of MOR, the spatial basis functions can be derived a 



priori using completeness conditions ( Ladyzhenskaia et al. 1969| Noack & Eckelmann 
1994[ ), from Navier-Stokes eigenfunctions ( Joseph||1976 ), or a posteriori from a snapshot 
of a solution data set, like the Proper Orthogonal Decomposition (POD) ( |Holmes et al. 
2012[ ) or Dynamic Mode Decomposition (DMD) ( |Rowley et aL|2009| |Schmid||2010[ ). 

The temporal coefficients, of the spectral discretized dynamical system are chosen 
such that the error is orthogonal to a subspace of H. Let {Vi € H\ i = 1, . . . ,n} be a 
basis for a subspace of H . We seek a.; such that 



- (vuX 



= 



(2.3) 



where (•, -) n is an appropriately defined spatial inner product. In the Galerkin projection 
approach, the trail basis are equal to the test basis; Vi = u t . The projection yields a set 
of evolution equations for the temporal coefficients 



fi(a) 



(2.4) 



Given a set of initial conditions the evolution equation can be integrated using standard 
numerical integration techniques. 



2.2. Proper orthogonal decomposition (POD) 

The Proper Orthogonal Decomposition (POD) is an information compression technique 
that is applicable to a wide variety of problems including digital image compres sion (|Richa -ds 
& Jia||2006| ), signal processing ( |Vaccaro||1991[ ) and bioinformatics fWall et a/.||2003[ ). We 
express the POD optimality condition in terms of the turbulent kinetic energy, e(t) of 
the flow defined as 

1 



e(f) 



|it(a;,t)| dx 



(2.5) 



POD provides spatial and temporal basis functions whose turbulent kinetic energy 

2 



(0 



^a{t)u(x) 



i=l 



dx 



(2.6) 



is as close as possible, on average, to the turbulent kinetic energy, e(t) of a solution 
snapshot of the Navier-Stokes, u s (x,t). In other words, the POD optimality condition 
can be formulated as follows 



arg mm 

a.u 



S.t. 



u s (x,t) - ^2ai(t)ui(x) 



dx 



(1) (u i ,u j ) n = 

(2) (ai,aj) T = 



/ Ui(x) ■ Uj(x)dx 



(2.7) 



T 



ai(t)aj(t)dt = Xi 



where Sij is the Kronecker delta and the diagonal values of the correlation matrix, A 



are the POD eigenvalues. The discretized equivalent of Eq. (2.7) is solved using the well 



known Eckart- Young theorem (Demmel 1997) of Singular Value Decomposition (SVD) 
and the method of snapshots ( Sirovich||1987 ). 
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2.3. Galerkin MOR of the Navier Stokes equation 

The evolution of the velocity field, u(x,t) of an incompressible, isentropic, Newtonian 
fluid is governed by the following Navier-Stokes equations 

V-it = (2.8a) 

— = vAu - V • (uu) - Vp (2.86) 



where v is the viscosity and p is the pressure. In the standard Galerkin MOR approach, 
the fluid velocity vector, u(x,t) is discretized spectrally using separable spatial, Ui(x) 
and temporal, di(t) basis functions 

n 

u(x, t) « M [0 '"' =UQ + ^2 a i(~t) U i( X ) ( 2 - 9 ) 

i=l 

where u is the temporal mean. Galerkin projection yields a set of n coupled, quadratic 
Ordinary Differential Equations (ODEs) in the following canonical form 

n n 

a%= ^ QijkdjCLk Lijdj + bi (2-10) 

j,k=l j=l 

For divergence-free spatial basis functions and steady dirichlet boundary conditions, the 
Galerkin matrices are of the form 

Qijk = (ui, V • (ujUk)) n (2.11a) 
Lij = {ui,vAuj) n + («<,-V • (u Uj)) Q + (iti,-V ■ (ujU )) n (2.116) 
bi = {ui,vAu Q - V • (it wo)) n (2.11c) 

For open flows the pressure term cannot be ignored. In practice however, the pressure 
term is rarely modeled in full. For example, the contribution from the pressure term can 



usually be modeled well using a linear fit of the linear Galerkin terms, (Galletti et al. 
2004) |Noack7FaI1[2005|. 



3. A Novel, Galerkin-based Model Order Reduction approach for the 
Navier-Stokes equations at high Reynolds number 

Turbulence is a flow regime characterized by chaotic, multi-scale dynamics in both 
space and time. POD basis functions are by construction, biased towards the large, 
energy containing scales of the turbulent flow. Galerkin ROMs generated using only the 
first most energetic POD basis functions are therefore not endowed with the natural 
energy dissipation of the smaller, lower energy turbulent scales. Since accounting for all 
the scales of the turbulent flow directly by including a large enough number of POD basis 
function is computationally prohibitive, the traditional approach as been to model the 
contribution of the small scales empirically. In this section, a new, non-empirical MOR 
approach for turbulent Navier-Stokes fluid flows is proposed. Instead of modeling the 
small, energy-dissipative scales using empirical eddy viscosity models, the small scales 
are modeled directly using spatial basis functions that are different from the standard 
POD basis functions. 

The productive/dissipative characteristics of a set of basis functions can be quantified 
using the time-averaged rate of turbulent kinetic production. The turbulent kinetic energy 
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of a spectrally discretized fluid flow is defined as 

f 



,[l,...,n] 



(*) 



y^ j a l {t)u l {x) 



1=1 



dx 



(3-1) 



Equation (3.1 ) can be expressed as a function of the temporal coefficients, ai(t) alone by 



exploiting the orthonormality of the spatial basis functions 

n 

e[1 ,..,»] (t) l£ ai(t) 2 



(3.2) 



A dynamical equation for the turbulent kinetic energy is derived as follows 



d e^(t) = ^5>(i) a = £«,(i)4i,(t) 



<7f 



i=l 



(3-3) 



Using Eq. (2.10) together with Eq. (3.3), we arrive at 



d n n n 

e [l,...,n]^_ ^ Q ijk a l a j a k + ^ l.,,u,ti, ■ ^"/i. 



(3.4) 



i,j',fe=l 



Finally, taking the time-average of Eq. (3.4), results in 



\dt 



(t)) — ^ Qijk^ijk+ ^ LjjXij+ bjHi 



(3.5) 



where T G I" xnx ", A <E E nx " and /i <E R" are matrices containing the temporal corre- 
lations 



Tijfe — (a,iaja,k) T 



(3.6a) 
(3.66) 
(3.6c) 



A positive rate of kinetic energy production, (^&- 1, '"' n '{t)) T > is associated with 
basis functions that resolve a greater proportion of large, energy-producing scales while 
a negative rate of kinetic energy production, (jjet lv "' n l(£)) < is associated with basis 
functions that resolve a greater proportion of small, energy-dissipative scales. When the 
basis functions are POD basis functions, the rate of kinetic energy production tends to be 
positive, i.e. (^e' 1 ' , "' n '(t)) T > 0. This is to be expected because the POD basis functions 
are biased toward the large, energy-producing scales of a turbulent flow. Galerkin-based 
ROMs of the Navier-Stokes equation deriving using only the first few most energetic POD 
basis functions tend to over predict the kinetic energy of the flow. The new proposed 
stabilization methodology achieves improved accuracy by deriving a new set of basis 
functions, a, Ui that model a greater proportion of small, energy-dissipation scales of 
the turbulent flow. The proposed stabilization methodology can be summarized as the 
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following minimization problem 



arg mm 
s.t. 



u s (x,t) -y^ai(t)ui(x) 



dx 



1) (Ui,Uj) n 

d 



(3.7) 



(2) 



dt 



e [l.-.n](f) 



where, now, the turbulent kinetic energy, e' 1, ■■■'"] (t) is evaluated using the new basis 
functions, a. Mi. The specific choice of the parameter e, referred to as the global dissi- 
pation parameter, is discussed in section 3.0.1 Since the objective function is infinitely 
dimensional, problem (3.7) appears to be intractable. However, the infinitely dimensional 
objective function can be reduced to a finite dimensional objective function by assuming 
the n new basis functions are spanned by the TV POD basis functions where N > n. In 
the spirit of related previous work ( Amsallem fc Farhat||2011 ), one can write 

(3.8a) 
(3.86) 



u = ux 

A = X T A 



where X £ M. Nxn is an orthonormal (X T X — I n xn) transformation matrix and the basis 
functions are vectorized and assembled in a matrix as follows 

(3.9a) 
(3.96) 
(3.9c) 
(3.9d) 



u = 


Ml 


u 2 ■ 


■■ u N ] 


u = 




U2 ■ 


■ ■ u n ] 




ai 


a 2 ■ ■ 


• a N \ 




ai 


a 2 ■ ■ 





Hence, the minimization problem (3.7) can be formulated as the following finite dimen- 
sional problem 



arg mm 
s.t. 



J \u s (x,t) - UXX T A\ 2 dx 



Hi 

{1)X T X = I nX n 



(3.10) 



(2) 



dt 



where the turbulent kinetic energy balance (i.e. constraint no. 2) is expressed as follows 



^ Qijk'iijk+ ^ LjjXjj+ bjfjj (3-11) 



ij,k=l 



The matrices Q £ R nxnxn , L £ R nxn and b £ K n are the Galerkin matrices associated 
with the new spatial basis functions, for i = 1, ■ ■ ■ , n, and T £ K nxrix ™ ; A £ R nxn 
and fl £ R n are matrices containing the new temporal correlations 



T ijk 



{aiaja k ) T 

(a l a j ) T 

(a,i) T 



(3.12a) 
(3.126) 
(3.12c) 
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These new matrices can be all expressed as a function of the transformation matrix, X 
as follows 



/v 



Q 



ijk 



XpiX^Qp;;X : k 

P =l 



i,j,k = !,-■■ ,n 



L = X l LX 



b = X T b 

N 

^ijk — J]] .\ ; „ A'.' T ; , :: A';/, i, j, k = 1, • 
P =i 

a = a: t aa: 
fi = x 1 V 



where Q e M. NxN * N , L g m nxtv and & g 
the POD spatial basis functions, Ui for i = 



(3.13a) 

(3.136) 
(3.13c) 

(3.13d) 

(3.13e) 
(3.13/) 

N are the Galerkin matrices associated with 
N. 



1 



To the best of the authors' knowledge, problem (3.10) does not admit any known 



closed form solutions and so, must be approximated numerically; details are provided 
in section 3.0.2 Furthermore, in its current formulation, solving problem (3.10) numer- 



ically would be computationally prohibitive and so further reductions in complexity are 
required. First, the objective function is reformulated. The objective function in prob- 



lem (3.10) measures the error between the snapshot solution, u s (x,t) and the spectral 
reconstruction, 'Y^ = idi{t)Ui{x) or, in matrix notation, UA = UXX T A. By definition, 
the reconstruction of the snapshot matrix using POD basis functions, J27=i a>i(t)iii(x) 
or, in matrix notation, UA is the optimal reconstruction. Where by optimal it is meant 
that no other combination of the temporal and spatial basis functions provide a greater 
resolution of time-averaged, turbulent kinetic energy of the snapshot matrix. In other 
words, the following inequalities always hold 



u s (x,t)\ dx ) ^ 

1? IT 



\UAfdx) ■ 

<! IT 



UA 



dx 



(3.14) 



Problem (3.10) can, therefore, be reformulated as follows 



N n 

argmin V X u - V (X T XX) i 
s.t. {l)X T X = I nxn 



(3.15) 



where 



\UA\*dx 
a It 



N 



UA 



dx 



E hi 

i=l 

n n 



(3.16a) 



(3.166) 



i=i 



i=i 



At this point, the second constraint of the minimization problem is reformulated. The 
most expensive calculation in the second constraint is the evaluation of the nonlinear 
Galerkin tensors, Q € jj nx ™ x ™ anc i Q g ^NxNxN ^Jq^ m odel the convective nonlinear- 
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ity of the Navier-Stokes equations. As in any fully spectral discretization, approximation 
and evaluation of nonlinearities can be computationally prohibitive. This is particularly 
true in the proposed approach since the tensor, Q must be evaluated at each iteration. 



Over the years, various approximations to tensor products have been derived (Carlberg 
et a/.|2011||Wang fc Ahuja|2008| ). In this paper, an alternative method to alleviate tensor 
computations is presented. It can be shown that for a large variety of flow bound ary con- 
ditions, the Galerkin tensor, Q is energy preserving, i.e. Ylijk Qijk { a iO-jO,k) T = (Kraich- 
nan fc Chen|1989 Cordier et aL||20i2} |Noack et aL||20lT ) . For flows where this condition 
holds, it can be shown that transformed Galerkin tensor Q is also energy preserving 
because a linear superposition of basis functions do not modify the boundary condi- 
tions. This nonlinear energy conversation principle holds for many boundary conditions 
including steady Dirichlct conditions, ambient flow or free-stream conditions at infinity, 
or periodic boundary conditions. For convective boundary condition, this energy con- 
servation principle has not been rigorously proven. Fortunately, it has been observed to 
hold in a variety of numerical simulations of flows with convective boundary conditions. 



Hence, for these boundary conditions, problem (3.10) can be reformulated as follows 



argmin V \ u - V (X T XX) 

xeR Nxr 



s.t. 



i=i i=i 
(l)X T X = I n 



(3.17) 



(2) J2 {X T LXoX T \X) ij =e 

where o is the Hadamard product. This final form of the minimization problem can be 
solved efficiently using standard numerical minimization algorithms. 

3.0.1. Global dissipation parameter, e 

For a statistically stationary turbulent flow, the rate of turbulent kinetic energy pro- 
duction is, on average, zero, i.e. (j^e) T = 0. It would seem plausible, therefore, that 
an accurate ROM could be derived using a set of basis functions, a and u for which 
(^e[ 1 ''"' rl l(t)) T = . Unfortunately, numerical simulations demonstrate that this not 
the case. Simulations suggest that ROM stability can be achieved using dissipative basis 
functions, i.e. basis functions for which (jjel 1 ''"'"] (i)) T = e where e is a small negative 
scalar. In this paper, ROM stability is measured using the statistics of the turbulent 
kinetic energy. Specifically, ROM stability is quantified using the first moment of the 
turbulent kinetic energy as follows 



'Mt? 



<5rom 



(3.18) 



where the temporal coefficients, 5j(t) in (3.18) are derived using numerical integration 
of the ROM and <5rom is the convergence criteria. Numerical simulations suggest that 
ROM stability is a smooth function of the dissipation parameter, e and therefore, a 
standard root-finding algorithm can be implemented. For both test cases summarized in 
this paper, convergence is achieved in less than 20 iterations that, depending on the size 
of the ROM, equates to several minutes of CPU time. 

3.0.2. The algorithm and numerical implementation 



For reasonably small problems, Eq. (3.17) can be directly implemented and solved 



in MATLAB using Sequential Quadratic Programming (SQP) (Dixon & Szego 1975 
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|Schittkowski|1986|[Nocedal fc Wright|1999[ ) and Brent's methods ( |Brent|2002| ; MATLAB 
functions fmincon and f zero respectively. Details of the MATLAB implementation are 
summarized in the appendix. The following pseudo-code summarizes the approach. 



Algorithm 1 Stabilization Algorithm 



Input: POD basis functions, Ui(x) and <Xj(i) for i = 1, ••• ,N and the associated 
Galerkin matrices Q e R NxNxN , L e R NxN and b G R N . The global dissipation 
parameter e, 

Output: Galerkin matrices Q € R" x " x ", L G R nxn and b e R n associated with the 
transformed basis functions, iii for i = 1, • • • , n 
1: Choose n > 0, s.t. N > n 
2: Compute the POD eigenvalues, Ay = (aiOj) T 



3: while ROM is unstable, i.e. 



Solve the optimization problem: 



> * 



ROM 



do 



/V 



argmin £ A. t - £ (A T AA% 
xeK«x„ , =1 , =1 

s.t. (1)A T A = I nxn 

n 

(2) ^ (A T iA oX T XX)..=e 



Evaluate new Galerkin matrices, Q, L and b. 
Integrate ROM, check for stability. 
Update e using Brent's method, 
end while 

Return u[, X and new Galerkin matrices Q, L and b 



4. Applications 

4.1. Lid-driven cavity 

The lid-driven cavity is a well known benchmark problem used to validate fluid flow nu- 
merical schemes and reduced order models ( Cazemier et q/.|1998 Shankar & Deshpande 



2000 Terragni et a/.|20TT ). Specifically, isothermal, incompressible, two-dimensional flow 



inside a square cavity driven by a prescribed lid velocity, una = (1 
The Reynolds number (Re) is defined with respect to the maximum velocity of the lid 
and the width of the cavity. The Navier-Stokes are discretized in space using Chebyshev 
polynomials. The convective nonlinearities are handled pseudo-spectrally and the Cheby- 
shev coefficients are derived using the Fast Fourier Transform (FFT) . The Navier-Stokes 
are discretized in time using a semi-implicit, second-order Euler scheme. Figure [l] is a 
snapshot of a statistically stationary solution at Re u = 3 x 10 4 . This particular simula- 
tion was performed using a 128 2 Chebyshev grid. The computations were performed on 
a cluster at Duke University using 8 processors. The simulation is first initialized over 
100 000 time steps (At = 1 x 10~ 3 s) which corresponds to a total run time of 1 hours. 
The data base is then generated: 50 000 iterations are done on 0.5 hours to create 5000 
snapshots. 



10 




11 


E [l,...,n 


1 


16.06 


2 


29.21 


3 


37.45 


4 


44.88 


5 


50.37 


10 


67.16 


20 


82.40 


50 


93.21 


200 


99.31 



Table 1. Percent of time-averaged, turbulent kinetic energy captured by the first n POD basis 
functions of the lid-driven cavity at Re u = 3 x 10 . 



4.1.1. POD basis functions of the lid-driven cavity 

A database of 5000 DNS snapshots of the lid-driven cavity were used to find the 
POD basis functions, itj(x) and a;(i). The nondimensional time interval between each 
snapshot equaled 0.1. Increasing the number of snapshots or the time interval between 
each snapshot had no significant effect on the performance characteristics of the ROMs. 
The percentaged of time-averaged, turbulent kinetic energy captured by the first n POD 
basis functions is labeled, ijl 1 -- ^"] and defined as follows 



E 



/ e [l.....n]\ 
[l,...,n] _ V s / 



^ x 100 = 



J2a(t)u(x) 



dx 



x 100 = 



\u s (x, t)\ 2 dx 



x 100 (4.1) 



Results for the lid-driven cavity are shown in Table [T] Vorticity contours of spatial POD 
basis functions, Ui(x) for i = 1, 2, 20, 50 and 200 of the lid-driven cavity are illustrated in 
Fig. [2] As expected, the overall spatial wavelengths of the POD basis functions tends to 
increases with order. The low-order POD basis functions correspond to the large, high- 
energy physical scales of the turbulent flows while the higher-order POD basis functions 
correspond to the small, low-energy physical scales of the turbulent flow. 
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Figure 2. Vorticity contours of spatial POD basis functions of the lid-driven cavity at 

Re u = 3 x 10 4 . 
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Figure 3. Galerkin ROMs of the lid-driven cavity derived using the first n = 5(£ 
10 < F — 1 >. 50( P^1 and 200(==| POD basis functions; DNS f^^j ) 



4.1.2. Galerkin ROMs of lid- driven cavity using standard POD basis functions, ai(t) 
and Ui(x) 

In this section, Galerkin ROMs of the lid-driven cavity are derived using the standard 
POD basis functions, ai{t) and Ui(x) for i = 1, • • • , n. The turbulent kinetic energy as 
predicted by these ROMs is illustrated in Figure [3] Numerical integration of the ROM 
was performed in MATLAB using built-in adaptive forth/fifth order Runge-Kutta scheme 
DDE45. Despite the fact that POD basis functions capture a large percentage of the time- 
averaged, turbulent kinetic energy of the snapshot solution (See Table[T]), Figure [3] clearly 
indicates Galerkin ROMs based on these basis functions do not converge to the DNS. 
Convergence is achieved only when a prohibitively large number of POD basis functions 
are utilized; approximately n = 200 POD basis functions for this particular test case. 
These convergence issues, of course, were anticipated. Galerkin ROMs of the Navier- 
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Figure 4. Time-averaged rate of turbulent kinetic energy production associated with the first 
n POD basis functions of the lid-driven cavity 



Stokes that utilize only the first few most energetic POD basis functions tend to under- 
resolve the small, energy-dissipating scales of the turbulent flow which, therefore, leads 
to excessive kinetic energy production. As discussed in Section [3j the time-averaged rate 
of kinetic energy production associated with the POD basis functions can be calculated 
using Eq. (3.5 1; these rates are illustrated in Figure |3j As shown, the time-averaged rate 



of kinetic energy production associated with the first n POD basis functions is always 
positive and approaches zero asymptotically. It is not surprising, therefore, that ROMs 
derived using only the most energetic basis functions tend to over predict the kinetic 
energy of the flow. 

4.1.3. Galerkin ROMs of lid- driven cavity using the new basis functions, di(t) and Ui(x) 

In this section, ROMs of the lid-driven cavity are derived using the new proposed 
methodology. Galerkin ROMs are derived using basis functions with negative rates of 
kinetic energy production, i.e. basis functions for which (^e[ 1 ''" ,rl l(t)) T = e where e is 
a negative scalar. In the previous section it was demonstrated that the standard POD 
basis functions have positive rates of kinetic energy production, (jje' 1 '"' ,n ''(i)) T > and 
therefore, produce Galerkin ROMs that over predict the kinetic energy of the flow. The 
critical rate of kinetic energy production, e for each ROM order n is found itcratively 



using the algorithm introduced in Section 3.0.2 



Figure [5] illustrates the evolution of the turbulent kinetic energy of lid-driven cavity 
as predicted by the new Galerkin ROMs. For the sake of brevity, only n = 5, 10 and 
20 ROMs are illustrated. The turbulent kinetic energy of the DNS solution is identified 
by the thick grey curves while the dashed and solid black curves identify the turbulent 
kinetic energy of the standard and new ROMs respectively. As predicted, ROMs derived 
using the new basis functions converge to the correct mean value of kinetic energy and the 
accuracy of the models is improved as n is increased. As stated previously, these new basis 
functions have negative rates of kinetic energy production, e while the standard POD 
basis functions have positive rates of kinetic energy production. Despite this difference, 
these new basis functions remain very similar to the POD basis functions as summarized 
in Table [2] The percentaged of time-averaged, turbulent kinetic energy captured by the 
first n new basis functions is labeled, E^ 1 *---'™] and defined as follows 

n ~ 

ft[l,...,n] = *=1 x 10Q ( 4 2 ) 

i=l 

For all orders n, the new basis functions capture a very similar total of kinetic energy 
of the flow. In other words, the transformation matrix associated with these new basis 
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n £[l-..n] 

5 50.37 50.01 
10 67.16 66.72 
20 82.40 82.11 

Table 2. Percentage of time-averaged, turbulent kinetic energy captured by the first n basis 

functions of the lid-driven cavity. 



functions resembles the following 



X 



Nxn 



{N - 



n) x n 



} Nxn 



(4.3) 



where /„ 



and O 



(N—n)xn 



are the identity and null matrices respectively. 5n xti is a 



matrix whose entries are all less than one, i.e. Sij < 1 Vi,j. Therefore, the new basis 
function inherent much of the optimality of the original POD basis function. 

4.2. Mixing layer 

The data base used for the present work corresponds to the DNS of an isothermal 2-D 
mixing layer. The numerical algorithm is the same as that employed previously for studies 
on jet noise sources ( Cavalieri et a/.|2011 ). The full Navier-Stokes equations for 2-D fluid 
motion are formulated in Cartesian coordinates and solved in conservative form. Spatial 
derivatives are computed with a fourth-order-accurate finite scheme for both the inviscid 



and viscous portion of the flux (Hayder & Turkel 1993 Gottlieb & Turkel 1976 ). A second- 



order predictor-corrector scheme is used to advance the solution in time. In addition, 
block decomposition and MPI parallelization arc implemented. The three-dimensional 
Navier-Stokes characteristic non-reflective boundary conditions (3D-NSCBC), developed 
in (Lodato et al. 20081, are applied at the boundaries of the computational domain to 



account for convective fluxes and pressure gradients across the boundary plane. In order 
to simulate anechoic boundary conditions, the mesh is stretched and a dissipative term is 
added to the equations, in the sponge zone ( Colonius et q/.|19"93 ). A detailed description 
of the numerical procedure is given in paviller||20f0| ~ 

The inflow mean streamwise velocity profile is given by a hyperbolic tangent profile 



u{y) = [/ 2 + -A[/[l+tanh(2j,)] 



(4.4) 



with AU = U\ — Ui the velocity difference across the mixing layer, where U\ and U2 are 
the initial velocity above and below, respectively. The velocities, lengths, and time are 
nondimensionalized with AU and the initial vorticity thickness S u . The flow Reynolds 
number is Re = 5 U AU /i> a = 500 where the subscript (•) indicates a constant ambient 
quantity. The Mach number of the free streams are Mi — U\jc a — 0.1 and M2 — C/2/ca = 
0.033, where c a is the sound speed. The inflow mean temperature is calculated with 
the Crocco-Busemann relation, and the inflow mean pressure is constant. The Prandtl 
number is selected to be Pr — 0.7. Finally, the convective Mach number is given by 
M c = AU/2c a = 0.033, so that the flow can be assumed as quasi- incompressible. The 
numerical code was extensively validated against numerical and experimental data; some 



results can be found in (Cavalieri et al. 2011 Daviller 2010). The computations were 



performed on the cluster of the PPRIME Institute in Poitiers, France using 64 processors. 
The computational domain comprises approximately 2.1 million grid points: 2367 points 
in the streamwise direction, and 884 points along the y direction. The extension of the 
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Figu re 5, E volution of the turbulent kinetic energy of the lid driven cavity as predicted by th e 
DNS {- - -\ , a standard POD-based Galerkin ROM ( | } and the new Galerkin ROM | |— [ ), 



computational domain is 325(5^ x 1205^. The sponge regions are from x = — 20<5 W to 
x = and x = 250(5^ to x = 305(5 W in the streamwise direction, and from ±5(W W to 
±604j in the transverse y direction. To promote a natural transition to turbulence from 
an initially laminar solution, the flow is forced by adding at every iteration solenoidal 
perturbations defined as ( |Bogey|2000 ).The simulation is first initialized over 330 000 time 
steps (At = 1.8 x 10 _8 s) which corresponds to a total run time of 35 hours. The data base 
is then generated: 1093 695 iterations are done on 115 hours to create 2000 snapshots. 



4.2.1. POD basis functions of the mixing layer 

A database of 2000 DNS snapshots of the mixing-layer cavity were used to find the 
POD basis functions, Ui(x) and a,i(t). The nondimensional time interval between each 
snapshot equaled 1. Increasing the number of snapshots or the time interval between each 
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Figure 6. Vorticity contours from a DNS of the mixing at Re^ = 500 



11 


E [l,...,n 


1 


18.05 


2 


34.69 


3 


40.20 


4 


4.62 


5 


50.71 


10 


66.80 


20 


81.77 


50 


93.93 


100 


98.36 



Table 3. Percent of time-averaged, turbulent kinetic energy captured by the first n POD basis 
functions of the mixing layer at Re^ = 500. 



snapshot had no significant effect on the performance characteristics of the ROMs. The 
percentaged of time-averaged, turbulent kinetic energy captured by the first n POD basis 
functions of the mixing layer are summarized in Table [T] Vorticity contours of spatial 
POD basis functions, Ui(x) for i = 1, 2, 20 and 100 of the mixing layer are illustrated in 
Fig. [7] As expected, the overall spatial wavelengths of the POD basis functions tends to 
increases with order. The low-order POD basis functions correspond to the large, high- 
energy physical scales of the turbulent flows while the higher-order POD basis functions 
correspond to the small, low-energy physical scales of the turbulent flow. 

4.2.2. Galerkin ROMs of mixing layer using standard POD basis functions, aj(t) and 

Ui(x) 

In this section, Galerkin ROMs of the mixing layer are derived using the standard POD 
basis functions, a,(£) and iii(x) for i = 1, • • • , n. The turbulent kinetic energy as predicted 
by these ROMs is illustrated in Figure [5J Similarly to the lid-driven cavity test case, 
Galerkin ROMs of the mixing layer based on POD basis function tend to over predict the 
kinetic energy of the flow. Similarly to the lid-driven cavity case, these inaccuracies were 
anticipated. As illustrated in Figure [9j the rate of kinetic energy production associated 
with the POD basis functions is positive. 

4.2.3. Galerkin ROMs of the mixing layer using the new basis functions, a^t) and Ui(x) 

In this section, the proposed new methodology is used to derive more accurate Galerkin 
ROMs of the mixing layer. Figure [10] illustrates the evolution of the turbulent kinetic 
energy of the mixing layer as predicted by the new Galerkin ROMs. For the sake of 
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Figure 7. Vorticity contours of spatial POD basis functions of the mixing layer at Re,^ = 500. 



brevity, only n = 5, 10 and 20 ROMs are illustrated. The turbulent kinetic energy of 
the DNS solution is identified by the thick grey curves while the dashed and solid black 
curves identify the turbulent kinetic energy of the standard and new ROMs respectively. 
As expected, the Galerkin-based ROMs of the lid-driven cavity derived using standard 
POD basis functions over predict the kinetic energy of the flow. On the other hand, 
ROMs derived using the new basis functions converge to the correct mean of kinetic 
energy of DNS. As stated previously, these new basis functions have negative rates of 
kinetic energy production, e while the standard POD basis functions have positive rates 
of kinetic energy production. Despite this difference, these new basis functions remain 
very similar to the POD basis functions as summarized in table [4] 
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Figure 9. Time-averaged rate of turbulent kinetic energy production associated with the first 
n POD basis functions of the mixing layer 



n E [1 '-' n] E [1 '-^ 
5 50.71 50.21 
10 66.80 66.67 
20 81.77 79.14 

Table 4. Percentage of time-averaged, turbulent kinetic energy captured by the first n basis 

functions of the mixing layer. 



5. Conclusions 

A new approach to Model Order Reduction of the Navier-Stokes Equations was pro- 
posed. The method provides spatial basis functions different from the usual proper 
orthogonal decomposition (POD) basis functions. These new basis functions resolve a 
greater range of physical scales of the turbulent fluid flow compared to the POD basis 
functions that are, by construction, biased toward the large, energy containing scales. 
When these basis functions are used in a Galerkin projection of the Navier-Stokes equa- 
tions, more stable and accurate ROMs are derived. 
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Appendix A. 

The following are MATLAB functions that solve the minimization problem for a given 
delta_e. The inputs of the function stabilizeROM are the N spatial basis functions, 
u and the temporal basis funcitons, a. The function stabilizeROM outputs the trans- 
formation matrix, X and matrices containing the new spatial basis functions, u_tilde 
and temporal basis functions a_tilde. The objective function, objective evaluates the 
difference between the optimal reconstruction of the time-averaged, turbulent kinetic en- 
ergy using the POD basis functions and the new basis functions. The constraint function 
constraint functions evaluates the time-averaged, turbulent kinetic energy. 



i function [X, u_t ilde , a_t i lde ] = stabilizeROM (u, a, L, N, n) 



2 global lambda L delta_e 

3 

4 for i=l:N; 

5 lambda ( i ) = mean ( a ( i , : ) . *a ( i, : ) ) ; 

6 end 

7 

8 xO = eye (n, n) ; 

g xO (N, : ) = 0; 

in 

n problem = createOptimProblem ( ' f mincon ' , . . . 

12 'objective', Sobjective, ... 

13 'nonlcon', ^constraint, ... 

14 ' xO ' , xO) ; 

15 [x, fval, EXITFLAG, OUTPUT, LAMBDA] = fmincon (problem) 
is OUTPUT. message 

17 

is X=x* (x ' *x) " ( — 1/2) ; 

19 u_tilde=u*X; 

2n a_tilde=X ' *a; 



2i end 



1 function [f] = objective (x) 

2 global lambda; 

3 X = x* (x ' *x) * (-1/2) ; 

4 sum_lambda = sum (lambda); 

5 lambda_tilda = diag (X ' *diag (lambda) *X) ; 

6 sum_lambda_tilde = sum ( lambda_t ilde ) ; 

7 

8 f = sum_lambda — sum_lambda_tilde; 

9 end 



1 


function [c,ceq] = constraint (x, delta_e ) 




2 


global D lambda; 




3 


X=x* (x' *x) " (-1/2) ; 




4 


L.tilde = X'*L*X; 




5 


lambda_tilde = (X ' *diag (lambda) *X) ; 




6 


ceq = sum (sum (D_tilde . *lambda_tilde) ) - 


- delta_e; 


7 


c = [ ] ; 




8 


end 
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